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The  aim  of  this  work  is  to  study  the  effects  of  gas-diffusion  layer  (GDL)  anisotropy  and  the  spatial  variation 
of  contact  resistance  between  GDLs  and  catalyst  layers  (CLs)  on  water  and  heat  transfer  in  polymer  elec¬ 
trolyte  fuel  cells  (PEFCs).  A  three-dimensional,  two-phase,  numerical  PEFC  model  is  employed  to  capture 
the  transport  phenomena  inside  the  cell.  The  model  is  applied  to  a  two-dimensional  cross-sectional  PEFC 
geometry  with  regard  to  the  in-plane  and  through-plane  directions.  A  parametric  study  is  carried  out 
to  explore  the  effects  of  key  parameters,  such  as  through-plane  and  in-plane  GDL  thermal  conductivi¬ 
ties,  operating  current  densities,  and  electronic  and  thermal  contact  resistances.  The  simulation  results 
clearly  demonstrate  that  GDL  anisotropy  and  the  spatial  variation  of  GDL/CL  contact  resistance  have  a 
strong  impact  on  thermal  and  two-phase  transport  characteristics  in  a  PEFC  by  significantly  altering  the 
temperature,  water  and  membrane  current  density  distributions,  as  well  as  overall  cell  performance.  This 
study  contributes  to  the  identification  of  optimum  water  and  thermal  management  strategies  of  a  PEFC 
based  on  realistic  anisotropic  GDL  and  contact-resistance  variation  inside  a  cell. 

©  2009  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

It  is  well  recognized  in  the  fuel  cell  community  that  water  man¬ 
agement  is  a  key  issue  in  the  operation  of  polymer  electrolyte  fuel 
cells  (PEFCs).  Without  adequate  balance  between  membrane  dehy¬ 
dration  and  water  flooding,  irreversible  losses  due  to  electrolyte 
or  mass  transport  resistance  become  excessive  and  exert  a  negative 
influence  on  the  performance  and  durability  of  PEFCs.  Furthermore, 
complex,  multi-dimensional,  two-phase  water  transport  phenom¬ 
ena  inside  PEFCs  are  strongly  affected  by  thermal  transport  due  to 
the  close  interaction  between  water  and  temperature,  where  a  sig¬ 
nificant  amount  of  heat  is  released  or  absorbed  due  to  the  phase 
change  of  water.  In  addition,  the  evaporation  and  condensation 
rates  of  water  are  also  a  strong  function  of  temperature.  Therefore,  it 
has  already  been  mentioned  in  the  literature  [  1  -3  ]  that  water  man¬ 
agement  in  PEFCs  should  be  coupled  with  thermal  management  to 
ensure  the  high  performance  and  durability  of  PEFCs. 

As  shown  in  the  scanning  electron  micrograph  of  Fig.  1 ,  a  typical 
carbon-paper  gas-diffusion  layer  (GDL)  in  PEFCs  is  comprised  of  car¬ 
bon  fibres  that  are  preferentially  oriented  in  the  in-plane  direction 
[4]  and  exhibit  significant  anisotropy  in  material  properties  along 
the  through-plane  and  in-plane  directions.  Therefore,  it  is  expected 
that  consideration  of  the  anisotropy  GDLs  in  a  fuel-cell  model  is  nec¬ 
essary  for  predicting  more  realistic  and  accurate  temperature  and 
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water  profiles  in  PEFCs.  To  date,  a  number  of  fuel-cell  models  of 
varying  levels  of  complexity  have  been  published  [  1  -3,5-20].  These 
range  from  simplified,  one-dimensional,  fuel-cell  component  mod¬ 
els  to  comprehensive,  computational  fuel-cell  dynamics  (CFCD) 
models  that  deal  with  multi-physics  and  multi-dimensional,  multi¬ 
phase  phenomena  in  fuel-cells.  By  contrast,  only  a  few  researchers 
[18-20]  have  attempted  to  investigate  the  effects  of  GDL  anisotropy 
on  water  and  heat  transport  phenomena.  Meng  [18]  introduced  a 
simplified  numerical  technique  to  approximate  anisotropic  elec¬ 
tron  transfer  phenomena  in  GDLs.  Once  the  in-plane  electronic 
conductivity  was  used  for  both  in-plane  and  through-plane  direc¬ 
tions  in  the  catalyst  layer  (CL)  and  GDL,  extra  contact  resistance  was 
applied  to  the  interface  between  a  GDL  together  with  a  bipolar  plate 
(BP)  to  compensate  for  the  effect  of  reduced  through-plane  elec¬ 
tronic  resistance.  Using  the  simplified  approach,  the  study  reported 
a  significant  improvement  in  numerical  convergence  and  stabil¬ 
ity  during  large-scale,  three-dimensional  simulations.  Pasaogullary 
et  al.  [19]  considered  anisotropic  heat  and  species  transport  in 
their  half-cell  model.  According  to  their  simulation  results,  the 
effect  of  GDL  anisotropy  is  indeed  significant  in  terms  of  the  con¬ 
siderable  change  in  the  water  and  temperature  distributions  in 
PEFCs.  Bapat  and  Thynell  [20]  presented  a  two-dimensional,  single¬ 
phase,  fuel-cell  model  and  investigated  the  effects  of  the  anisotropic 
thermal  conductivity  of  a  GDL  and  membrane,  and  the  thermal 
contact  resistance  at  the  interface  between  the  GDL  and  the  bipo¬ 
lar  plate.  They  reported  the  considerable  effects  of  the  anisotropic 
thermal  conductivity  of  a  GDL  on  the  temperature  distributions  in 
PEFCs  and  concluded  that  a  high,  in-plane  thermal  conductivity  is 
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Nomenclature 

a  water  activity  or  effective  catalyst  area  per  unit  of 

total  volume  (m2  m-3) 

A  area(m2) 

Cp  specific  heat  (kj  kg-1 I<-1 ) 

C  molar  concentration  (mol  m-3 ) 

Dk  mass  diffusivity  of  species  k  (m2  s-1 ) 

EW  equivalent  weight  of  dry  membrane  (kg  mol-1 ) 

F  Faraday  constant  (96,487  C  mol-1 ) 

i0  exchange  current  density  (A  m-2 ) 

j  transfer  current  density  (A  m-3 ) 

jl  diffusive  mass  flux  of  i  phase  (kg  m-2  s-1 ) 

J  Leverett  function 

h  enthalpy  per  unit  mass  (kj  kg-1 ) 

k  thermal  conductivity  (W  m  K-1 ) 

kr  relative  permeability 

I<  hydraulic  permeability  (m2 ) 

M  molecular  weight  (kg  mol-1 ) 

m  mass  fraction 

n  number  of  electrons  in  electrochemical  reaction  or 

diffusivity  correction  factor 
nc  catalyst  coverage  coefficient 

nd  electroosmotic  drag  coefficient 

P  pressure  (Pa) 

Pc  capillary  pressure  (Pa) 

RH  relative  humidification  of  the  inlet 
Ru  universal  gas  constant  (8.314J  mol-1 I<-1 ) 
s  stoichiometry  coefficient  in  electrochemical  reac¬ 

tion  or  liquid  saturation 
S  source  term  in  transport  equation 

t  time  (s) 

T  temperature  (I<) 

u  fluid  velocity  and  superficial  velocity  in  porous 

medium  (ms-1) 

U0  thermodynamic  equilibrium  potential  (V) 

V  volume  (m3) 

Vcen  cell  potential  (V) 

Greek  letters 

a  transfer  coefficient 

8t  thickness  of  component  i 

s  volume  fraction  of  gaseous  phase  in  porous  region 

se  volume  fraction  of  ionomer  phase  in  CL 

y  advection  correction  factor 

X  membrane  water  content  ( mol  H2  0 / mol  S03  - ) 

Xa  relative  mobility  of  phase  a 

0  phase  potential  (V) 

1 7  overpotential  (V) 

0  contact  angle  (°) 

fi  viscosity  (kg  m-1  s-1 ) 

p  density  (kg  m-3) 

pmem  membrane  density  ( kg  m-3 ) 

v  kinematic  viscosity  (m2  s-1 ) 

a  surface  tension  (Nm-1)  or  electronic  conductivity 

(Sm-1) 

r  viscous  shear  stress  (N  m-2) 

k  ionic  conductivity  (S  m-1 ) 

§  stoichiometry  flow  ratio 

Superscripts 
e  electrolyte 

eff  effective  value  in  porous  region 

mem  membrane 


8  gas 

l  liquid 

ref  reference  value 

s  solid 

sat  saturation  value 

Subscripts 
a  anode 

avg  average  value 

BP  bipolar  plate 

c  cathode 

CL  catalyst  layer 

e  electrolyte 

g  gas  phase 

GC  gas  channel 

GDL  gas-diffusion  layer 

H2  hydrogen 

i  species  index 

in  channel  inlet 

m  mass  equation 

mem  membrane 

N2  nitrogen 

02  oxygen 

ref  reference  value 

t  total 

s  solid 

sat  saturation  value 

w  water 

0  potential  equation 

0  standard  condition  (298.15  K  and  101.3  kPa  (1  atm)) 


desirable  as  it  leads  to  more  uniform  temperature  distribution 
inside  PEFCs. 

Although  these  previous  studies  [18-20]  demonstrated  the  sig¬ 
nificant  influences  of  GDL  anisotropy,  some  profound  effects  in  an 
anisotropic  GDL  have  not  been  studied  in  detail.  These  include  the 
inter-relationship  between  thermal  gradients  and  two-phase  water 
transport,  and  the  two-phase  heat  transfer  that  is  driven  by  the 
heat  release/absorption  that  arises  from  the  phase  change.  In  the 
present  work,  anisotropic,  two-phase  heat  and  water  transport  phe¬ 
nomena  in  PEFCs  are  investigated  by  implementing  the  anisotropic 
functions  into  a  three-dimensional  two-phase  PEFC  model  that  has 
been  developed  earlier  [16,17].  In  addition  to  studying  anisotropic 
heat  and  water  transport,  an  investigation  is  made  of  the  effects 
of  thermal  contact  resistance  at  both  the  GDL/CL  and  the  GDL/BP. 
Previously,  the  effect  of  thermal  contact  resistance  was  studied  by 
Bapat  and  Thynell  [20],  who  applied  a  constant  value  to  the  contact 


Fig.  1.  SEM  image  of  TGP-H  carbon-paper  GDL  [4]. 
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resistance  at  the  interface  between  a  GDL  and  BP  but  neglected  the 
contact  resistance  at  the  interface  between  a  GDL  and  CL.  In  reality, 
the  contact  resistance  between  a  GDL  and  CL  is  expected  to  be  more 
significant  and  vary  along  the  in-plane  direction  between  the  gas 
channel  (GC)  and  land  regions.  This  is  because  the  assembly  clamp¬ 
ing  pressure,  i.e.,  the  compression  pressure  that  results  from  torque 
on  the  bolts  that  clamp  the  fuel  cell,  cannot  be  fully  applied  to  GDL 
areas  under  the  gas  channels.  Therefore,  in  this  study,  the  effect 
of  the  variation  of  the  contact  resistance  at  the  interface  between 
a  GDL  and  CL  along  the  in-plane  direction  is  considered  for  pre¬ 
dicting  more  realistic  water  and  temperature  profiles  inside  a  PEFC. 
This  study  assists  in  finding  appropriate  water  and  thermal  man¬ 
agement  strategies  of  a  PEFC  under  realistic  anisotropic  GDL  and 
contact-resistance  variation. 

2.  Numerical  model 


Charge  conservation 

proton  transport :  V  •  (/ce^V0e)  +  5^  =  0 

(6) 

electron  transport :  V  •  (ae^V0s)  -S&  =  0 

(7) 

Energy  conservation 

V  •  (puCpT)  =  V  ■  (k^VT)  +  V  •  ( h°gplul )  +  St 

(8) 

In  the  conservation  equations,  Eqs.  ((1  )-(8)),  the  two-phase  mix¬ 
ture  properties  are  defined  as  a  function  of  s  and  (1  -s),  which 
denote  the  volume  fraction  of  the  open  pore  spaces  that  are  occu¬ 
pied  by  the  liquid  and  gaseous  phases,  respectively  [21  ]. 

Density  [21] 

p  =  pl  ■s  +  pg  ~s)  (9) 


2.1  Multi-dimensional  fuel-cell  model 

The  two-phase,  non-isothermal  PEFC  model  used  for  this  study 
is  based  on  the  multi-phase  mixture  (M2)  model  developed  by 
Wang  and  Cheng  [21].  The  fuel-cell  model  has  been  described  in 
detail  in  previous  publications  [16,17];  thus,  only  a  brief  summary 
of  the  conservation  equations  is  repeated  here,  as  follows. 

2.2.  Model  assumptions 


In  Eq.  (9),  the  gas  mixture  density,  pg ,  which  is  described  by  the 
ideal  gas  law,  varies  with  the  composition  (denoted  here  by  the 
mass  fractions,  {m^}).  That  is, 

^  =  S>fM  (10) 

Velocity  [21] 

pu  =  plul  +  pgug  (11) 


Employing  the  M2  formulation  to  describe  two-phase  transport 
in  PEFCs,  the  present  two-phase  PEFC  model  makes  the  following 
assumptions:  (1)  ideal  gas  mixtures  in  a  single-phase  region;  (2) 
laminar  flow  due  to  small  flow  velocities;  and  (3)  two-phase  mist 
flow  (i.e.,  homogeneous  flow)  in  gas  channels,  assuming  that  chan¬ 
nel  flooding  is  minimal  where  tiny  droplets  exist  and  travel  with 
the  gas  flow  inside  the  gas  channels. 

2.3.  Conservation  equations 


Relative  permeability 
k'  =  s4 


fc?  =  (l-s)4 

Kinematic  viscosity  [21  ] 


'k[  kf 

v=[f  +  ^ 


-1 


(12a) 

(12b) 


(13) 


Under  the  above  assumptions,  the  two-phase  PEFC  model  is  gov¬ 
erned  by  the  conservation  of  mass,  momentum,  chemical  species, 
charge,  and  thermal  energy. 

Mass  conservation 

V  •  (pu)  =  Sm  (1) 

Momentum  conservation 

flow  channels  (Navier-Stokes  Eqs.) 

-L- V  •  (pm)  =  -Vp  +  V  •  r  (2) 

£z 

porous  media  (Darcy’s  Eq.) 

pu  =  -*Vp  (3) 

Species  conservation 
flow  channels  and  porous  media 

V  •  ( YiprriiU )  =  V  .  [p«Dfe//V(mf )]  +  V  •  [(mf  -  m']/']  +  S,  (4) 


In  Eq.  (13),  Vs  is  the  kinematic  viscosity  of  a  gas  mixture  that 
varies  with  the  gas  composition  [22] 


^4  =  ^ 


XjPi 


pS  pS^Y"  XQ 

1=1  J  v 


,  where 


*  1  (,  Mi 

<Z>„  =  -=  1  +  tt- 


V8 


-1/2 


-44 


1/4' 


and 


/Zj[Nsm-2]  = 


'  ph2  =  o.2i  x  io-6r0-66 

fiw  =  0.00584  x  10-6T129 
<  pN.2  =  0.237  x  10-6T°-76 
„  p02  =  0.246  x  10-6T°-78 
T  in  Kelvin 


(14) 


Relative  mobility  [21  ] 


v.l, 


(15a) 


water  transport  in  membrane 


A*  =  1  -  V 


(15b) 


(nmem  \  /  j  \  /  ismem  \ 

WD™mV4  Mw  -  V  •  (nd-)  Mw  +  V  .  (—VP')  =  0 

(5) 


The  species  diffusivity  in  the  gas  mixture,  Df,  in  the  first  term 
on  the  right-hand  side  of  Eq.  (4)  is  defined  as  follows,  in  order  that 
the  summation  of  interspecies  diffusion  within  the  gaseous  phase 
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is  equal  to  zero  [22]. 

D8  =  — — — — - ,  where 

vj  ;  xj/Dfj 
j  +  i 

1.01310  ~7.p-75  f  1  1  V/2  (16^ 

iJ=  p  (x,i/3+x;/3)2  \^+^J 

Xhw  =  7.07,  Xw  =  12.7,  Xnw  =  17.9,  Xow  =  16.6 

Note  that  the  gaseous  diffusive  transport  for  a  porous  medium 
can  be  controlled  by  the  Knudsen  diffusion  effect  due  to  molecular- 
to-wall  collision  as  well  as  molecular  diffusion  that  is  caused 
by  molecular-to-molecular  collision,  as  described  in  Eq.  (16).  The 
Knudsen  diffusion  coefficient  can  be  computed  according  to  the 


The  transport  properties  of  electrolytes  are  correlated  with  the 
water  content  of  the  membrane,  X,  which,  in  turn,  is  a  function  of 
the  water  activity,  a,  as  follows  [24]. 


CjRuT 

Psat 


(23) 


X= 


Xs =0.043  +  17.81a  -  39.85a2  +  36.0a3  forO  <  a  <  1 
X'  =  22 


(24) 


The  electro-osmotic  drag  coefficient,  nd,  the  proton  conductiv¬ 
ity  in  the  membrane,  /c,  and  the  water  diffusion  coefficient  in  the 
membrane,  D™em,  have  been  reported  by  Springer  et  al.  [24]. 


nd  = 


2.5X 

22“ 


k  =  (0.51391  -0.326)  exp 


1268 


V303 


i)] 


(25) 

(26) 


nmem 

uw 


'  2.692661843  ■  10“10 

{0.87(3  -  X)  +  2.95(X  -  2))  ■  10"10  •  e(7-9728-(2416/T)) 

{2.95(4  -  X)  +  1 ,642454(X  -  3))  •  10"10  ■  e(7-9728-(2416/7)) 

„  (2.563  -  0.33X  +  0.0264X2  -  0.000671X3)  •  10“10  ■  e(7-9728-(2416/7')) 


fork  <  2 
for  2  <  X  <  3 
for 3  <  X  <4 
for  4  <  X  <  Xga=1 


(27) 


kinetic  theory  of  gases  as  follows. 


Df  = 


2  /8RUT\1/2 

3  V  JTMi  )  Yp 


(17) 


The  effective  diffusivity  of  species  i  in  the  gas  mixture  is  then 
obtained  by  combining  both  molecular  and  Knudsen  diffusion 
effects  with  the  effects  of  porosity  and  tortuosity  of  the  porous 
medium  by  using  the  Bruggeman  correlation  [23]. 


Df’e-ff  =  [s(l-s)]n 


1  1 

D?  +  Df 


-1 


(18) 


In  Eq.  (18),  it  is  seen  that  the  effective  gas  diffusivity,  D fe^,  in 
the  two-phase  region  is  a  function  of  both  the  porosity,  £,  and  the 
liquid  saturation,  s. 

On  the  other  hand,  the  diffusive  mass  flux  of  the  liquid  phase, 
j1,  which  is  shown  in  the  second  term  on  the  right-hand  side  of  Eq. 
(4),  is  expressed  as  a  function  of  the  capillary  pressure,  Pc,  such  that 
[21]: 


jl  =  plul  -  Xlpu  =  ^ XlXgVPc 


(19) 


2.4.  Source  terms  and  physicochemical  relations 

In  the  conservation  equations,  Eqs.  (( 1  )-(8)),  Sm,  Sit  S&  and  ST  that 
are  summarized  in  Table  1  denote  their  corresponding  source/sink 
terms  for  mass,  species,  charge  and  heat,  respectively.  These  source 
terms  result  from  the  electrochemical  reactions,  the  hydrogen  oxi¬ 
dation  reaction  (HOR)  in  the  anode  CL,  and  the  oxygen  reduction 
reaction  (ORR)  in  the  cathode  CL,  respectively,  which  are  repre¬ 
sented  by  kinetic  expressions  as  follows. 

/  r  \  1/2 

HOR  in  the  anode  CL  :  j=(l  -  s)ncairfa  (  J  ( 

(36) 


ORR  in  the  cathode  CL  : 


j  =  -(1  -s)  cai( 


fcniref 
o,c 


Cq2 


-02,ref 


3/4 


x  exp 


(37) 


In  Eq.  (19),  the  correlation  of  the  capillary  pressure,  Pc,  and  the 
Leverett  function  J(s),  is  given  by: 

,  f£\^2 

Pc=Pg  -P1  =cr cos0(—j  J(s)  (20) 

J  1.417(1 -s)- 2.120(1 -s)2  +  1.263(1  -s)3  if0c  <  90° 

\  1.417s -2.120s2  +  1.263s3  if0c  >  90° 

(21) 


The  second  term  on  the  left-hand  side  of  the  species  equation 
represents  the  advective  term,  in  which  the  advection  correction 
factor,  yi,  is  given  by  Wang  and  Cheng  [21  ],  i.e., 


p(Xlm\  +  X^mf) 

( splm\  +  (1  -  s)psm :?) 


(22) 


In  Eqs.  (36)-(37),  (1  -  s)nc  approximates  the  effect  of  decreasing 
electrochemically  active  catalyst  sites  due  to  the  presence  of  water 
in  the  CLs. 

The  surface  overpotential  in  Eqs.  (36)-(37)  is  defined  by 

HOR  in  the  anode  CL  :  rj  =  0S  -  4>e  (38) 

ORR  in  the  cathode  CL  :  p  =  0S  -  0e  -  H0  (39) 

In  Eqs.  (38)-(39),  the  thermodynamic  equilibrium  potential 
between  the  cathode  and  anode  is  given  by  [25] 

Uo  =  1.23-  0.9  x10-3(T-  298.15)  (40) 


Based  on  the  experimental  data  of  Parthasarathy  et  al.  [26],  the 
temperature  dependence  of  the  ORR  kinetic  parameter  is  approxi¬ 
mated  as  follows: 


af0efc(T)  =  air0efc( 353  K)-exp 


Ea 

P  1  \1 

.  Ru 

\T  353. 15  A 

(41) 


Table  1 

Two-phase  steady-state  PEFC  model:  source/sink  terms. 
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Source/sink  terms 


Mass 

In  the  CLs:  Sm  =  Y\ St+M WV  ■  (D™em(pmem /EW)VX) 

(28) 

Species 

For  water  in  CLs:  S,  =  M,[- 

-V  •  ((iid/f)/)  -  (sJ/nF)] 

(29a) 

For  other  species  in  CLs:  S 

h  =  -M,(sj/nF) 

(29b) 

Charge 

In  the  CLs:  S<&  =j 

(30) 

Energy 

In  the  CLs:  ST  =j(r]  +  T(dU0/dT ))  +  (BJk^) 

(31) 

In  the  membrane:  Sj  =  —e 
1 

ff 

f  Mi  =  chemical  formula  of  species  i 

(32) 

Electrochemical  reactions 

\  s,-M?  =  ne~ ,  where  < 
,  '  1 

Si  =  stoichiometry  coefficient 
[  n  =  number  of  electrons  transferred 

(33) 

Hydrogen  oxidation  reaction  (HOR)  in  anode  side:  H2-2H+  =2e~ 

(34) 

Oxygen  reduction  reaction  (ORR)  in  cathode  side:  2H20  -  02  -  4H+=4e~ 

(35) 

2.5.  Boundary  conditions  and  numerical  implementation 

The  inlet  velocities  in  the  anode  and  cathode  gas  channels  can 
be  expressed  as  a  function  of  the  anode  and  cathode  stoichiometric 
ratios  (§fl  and  respectively),  the  PEFC  operating  current  density 
(/),  the  cross-sectional  areas  of  the  anode  and  cathode  gas  chan¬ 
nels  ( Aa  and  Ac,  respectively),  and  the  densities  of  hydrogen  or  air 
gas  which  are  respectively  a  function  of  the  anode/cathode  inlet 
pressure,  temperature  and  humidity. 

tja{I/2F)AmemMn2  §c(f/4F)i4memM  q 

Win, a  = - T— S - 1  and  uinx  = - — — - 1  (42) 

PH2Aa  Po2Ac 

In  addition,  the  isothermal  boundary  condition  (80  °C  for  this 
study)  is  applied  to  the  outer  surfaces  of  the  bipolar  plate  for  tem¬ 
perature  calculations. 

The  two-phase  PEFC  model  described  above  is  numerically 
implemented  with  a  commercially  available  computational  fluid 
dynamics  (CFD)  package,  i.e.,  STAR-CD,  based  on  its  user-subroutine 
files  [27]. 

3.  Results  and  discussion 

3.1  Simulation  cases 

In  this  study,  the  comprehensive  two-phase,  non-isothermal 
PEFC  model  is  applied  to  two-dimensional  cross-section  geome¬ 
try  as  displayed  with  its  mesh  configuration  in  Fig.  2,  where  fully 
humidified  (100%)  hydrogen  and  air  gases  are  assumed  along  the 
anode  and  cathode  GCs.  The  properties  of  the  individual  cell  com¬ 
ponents  and  operating  conditions  are  listed  in  Tables  2  and  3, 
respectively,  and  the  relevant  physical  properties  are  summarized 
in  Table  4.  Using  the  cell  configuration,  properties  and  operating 
conditions,  a  parametric  study  has  been  conducted  to  investigate 
the  effects  of  GDL  anisotropy  and  contact  resistance  on  the  liq¬ 
uid  saturation  and  temperature  distribution,  as  well  as  the  overall 


Table  3 

Operating  conditions. 


Description 

Value 

Anode  inlet  pressure  (Pa) 

1.5  atm 

Anode  stoichiometry  (fa) 

3 

Cathode  inlet  pressure  (Pc) 

1.5  atm 

Cathode  stoichiometry  (£c) 

2.3 

RH  of  anode  inlet  ( RHa ) 

100% 

RH  of  anode  inlet  (RHC) 

100% 

cell  performance.  Based  on  typical  thermal-conductivity  values  of 
anisotropic  GDLs  [4,19,20],  the  thermal  conductivity  of  GDLs  is 
assumed  to  be  10.0  W  nrr 1  K-1  along  the  in-plane  direction  and  the 
through-plane  thermal-conductivity  value  is  taken  to  be  an  order 
of  magnitude  smaller  than  the  in-plane  value  (thus,  1.0  Wm-1  K-1 
along  the  through-plane  direction). 

In  terms  of  the  GDL  thermal  conductivities,  three  types  of  GDL 
are  defined  for  this  study.  The  first  GDL  is  tailored  to  have  a 
through-plane  thermal  conductivity  of  1.0  Wm-1  K-1  in  all  direc¬ 
tions  (denoted  the  ‘thru-plane  GDL’  in  the  following  discussion).  For 
the  second  GDL,  the  in-plane  value  of  10.0 Wm-1  K-1  is  assumed 
for  all  directions  (denoted  the  ‘in-plane  GDL’).  The  third  GDL  is 
anisotropic  and  has  different  thermal  conductivity  values,  viz., 
1.0  Wm-1  K-1  for  the  through-plane  direction  and  10.0  Wm-1  K-1 
for  the  in-plane  direction  (denoted  the  ‘anisotropic  GDL’). 

The  effects  of  the  contact  resistance  on  heat  and  electron  trans¬ 
port  inside  PEFCs  has  also  been  investigate.  As  shown  in  Table  2, 
constant  values  of  electronic  and  thermal  contact  resistance  are 
applied  to  the  GDL/BP  interfaces.  On  the  other  hand,  as  mentioned 
earlier,  the  spatial  variation  of  the  thermal  and  electronic  contact 


Table  4 

Physical  properties. 


Description 

Value 

Exchange  current  density  x  ratio  of  reaction  surface  to 

1.0  x  109  Am-3 

CL  volume  in  anode  side  (ai[J^) 

Exchange  current  density  x  ratio  of  reaction  surface  to 

1.0  x  104  Am-3 

CL  volume  in  cathode  side  (aij^) 

Activation  energy  for  oxygen  reduction  reaction  in 

73,269  J  mol-1 

cathode  side  (Ea) 

Reference  hydrogen  molar  concentration  (cH2,re/) 

40.88  mol  m-3 

Reference  oxygen  molar  concentration  (co2,re/) 

40.88  mol  m-3 

Anodic  and  cathodic  transfer  coefficients  for  hydrogen 

aa=oic  =  1 

oxidation  reaction  (HOR) 

Cathodic  transfer  coefficient  for  oxygen  reduction 

ofc  =  l 

reaction  (ORR) 

Dry  membrane  density  ( pmem ) 

2000  kg  m-3 

Equivalent  weight  of  electrolyte  in  membrane  ( EW ) 

1.1  kg  mol-3 

Faraday  constant  (F) 

96,487  C  mol”1 

Universal  gas  constant  (Ru) 

8.314  J  moP1  K-1 

Surface  tension  (or) 

0.0625  Nm-1 

Liquid  water  density  ( pl  (80  °C)) 

972  kg  irr3 

Liquid  water  viscosity  (pi1) 

3.5  x  10-4  N  s  m-2 

Table  2 

Cell  properties. 


Description  Value 


Porosity  of  GDLs  (£gdl) 

0.6 

Porosity  of  CLs  (£Cl) 

0.6 

Volume  fraction  of  ionomer  in  CLs  (£mc) 

0.18 

Permeability  of  GDLs  (I<gdl) 

3.0  x  10-12 

m2 

Hydraulic  permeability  of  membrane  (Kmem) 

5.0  x  10-20 

m2 

Contact  angle  of  GDLs  and  CLs  (9) 

110° 

Effective  electronic  conductivity  in  CL  (a cl) 

1000  Sm”1 

Effective  electronic  conductivity  in  GDL  (crGDL) 

10,000  Sm 

-i 

Effective  electronic  conductivity  in  BP  (<xBP) 

20,000  Sm 

-i 

Electric  contact  resistance  between  BP  and  GDL  (Re,Bp) 

5.3  mQ,  cm; 

i 

Thermal  conductivity  of  BP  ( ksp ) 

20.0  Wm-1 

K-1 

Thermal  conductivity  of  membrane  (kmem) 

0.95  Wm”1 

K-1 

Thermal  contact  resistance  between  BP  and  GDL  (Rt,Bp) 

0.75  x  10“4 

cm2 
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are  applied 


Cathode  inlet 

Fig.  2.  Dimension  of  the  PEFC  two-dimensional  cross-section  geometry  (a)  and  its  mesh  configuration  (b)  where  the  thermal  and  electronic  contact  resistances  are  applied 
at  the  GDL/CL  interfaces  and  BP/GDL  interfaces. 


resistance  at  the  GDL/CL  interfaces  is  considered  along  the  in-plane 
direction  due  to  the  differential  clamping  pressures  that  are  applied 
at  the  gas  channel  and  land  regions,  respectively.  According  to 
the  experimental  contact-resistance-pressure  constitutive  relation 
provided  by  Zhang  et  al.  [28],  the  contact  resistance  between  the 
BP  and  the  GDL  exponentially  increases  with  decreasing  champing 
pressure.  Based  on  the  experimental  results,  it  is  concluded  that 
the  contact  resistance  at  the  GDL/CL  interface  under  the  channel 
region  is  exponentially  related  to  the  distance  from  the  edge  of  the 
land/channel  as  follows: 

Rchan=Rlande3-5x  (43) 

In  Eq.  (43),  x  denotes  the  distance  from  the  land/channel  inter¬ 
face  towards  the  channel  region  along  the  in-plane  direction  and 
the  factor  3.5  is  arbitrarily  assumed.  Based  on  the  above  assump¬ 
tion,  the  electronic  and  thermal  contact  resistances  at  the  GDL/CL 
interface  are  plotted  in  Fig.  3  along  the  in-plane  direction. 

3.2.  Simulation  results 

Fig.  4  displays  the  current  density  (/)  and  temperature  (T)  pro¬ 
files  in  the  membrane  along  the  in-plane  direction  for  three  types 
of  GDL,  where  PEFC  operation  at  1.5  A  cm-2  is  considered  and  the 
effects  of  both  electronic  and  thermal  contact  resistances  at  the 
GDL/CL  interfaces  are  ignored.  Since  the  heat  that  is  generated  in 
PEFC  is  primarily  removed  through  the  GDL  to  the  current-collector 
land  of  BP  by  lateral  conduction  (and  thus,  the  land  acts  as  a  heat¬ 
sink  for  the  GDL),  the  main  heat-flow  path  is  from  inside  CLs  toward 
land  regions.  Consequently,  for  all  three  cases,  higher  temperature 
is  predicted  near  the  channel  region  due  to  a  lack  of  effective  cool¬ 


ing  in  the  channel  region.  While  the  thru-plane  GDL  case  leads 
to  the  greatest  rise  in  membrane  temperature  due  to  its  low  GDL 
thermal  conductivity  (l.OWm-1  K-1),  the  difference  between  the 
temperature  profiles  of  the  in-plane  and  anisotropic  GDL  cases  is 
relatively  small.  This  indicates  that  an  isotropic  treatment  using 
the  in-plane  GDL  thermal  conductivity  leads  to  a  smaller  deviation 
in  the  temperature  profile  than  the  use  of  the  thru-plane  value. 
While  the  rise  in  the  membrane  temperature  that  corresponds  to 
the  anisotropic  GDL  is  calculated  to  be  roughly  2.1  °C,  a  rise  of  4.1  °C 
in  the  maximum  membrane  temperature  is  predicted  for  the  thru- 
plane  GDL,  which  demonstrates  a  significant  over-prediction  of  the 
rise  in  the  membrane  temperature  in  that  case.  On  the  other  hand, 
the  membrane  temperature  rises  up  to  1.3  °C  for  the  in-plane  GDL, 
which  implies  that  the  degree  of  under-prediction  of  the  rise  in  the 
membrane  temperature  by  the  use  of  the  in-plane  thermal  conduc¬ 
tivity  is  relatively  less  serious.  In  addition,  it  is  seen  in  Fig.  4  that 
the  shape  of  the  current  density  distribution  in  the  membrane  is 
considerably  altered,  when  the  rise  in  the  local  membrane  temper¬ 
ature  is  significant.  For  both  the  in-plane  and  anisotropic  GDLs,  the 
current  density  is  greater  near  the  channel  region.  This  indicates 
that  the  current  density  for  both  cases  is  mainly  determined  by 
the  limitation  in  oxygen  transport.  In  contrast  to  these  two  cases, 
the  current  density  near  the  channels  is  lower  for  the  thru-plane 
GDL  and  implies  that  significant  membrane  dehydration  occurs  due 
to  the  considerable  over-prediction  of  the  membrane  temperature 
near  the  gas  channels.  Thus,  ohmic  loss  is  a  dominant  factor  in  the 
membrane  current  density  and  cell  performance  for  the  thru-plane 
GDL  case. 

Fig.  5  presents  the  current  density  and  membrane  temperature 
profiles  for  the  same  three  GDL  cases  at  1.5  A  cm-2  but  the  spatial 
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Fig.  3.  Thermal  contact  resistance  (a)  and  electric  contact  resistance  (b)  variations 
at  the  GDL/CL  interfaces  along  the  in-plane  direction. 


variation  of  electronic  and  thermal  contact  resistances  between  the 
GDL  and  CL  is  taken  into  account.  A  comparison  of  Fig.  5  with  Fig.  4 
clearly  shows  that  the  presence  of  contact  resistance  between  the 
GDL  and  the  CL  (plotted  in  Fig.  3),  results  in  a  significant  rise  in 
temperature,  particularly  near  the  channel  regions.  For  each  GDL, 
the  rise  in  temperature  is  predicted  to  be  roughly  2°C  higher  in 
Fig.  5  than  in  Fig.  4.  The  thermal  effect  also  significantly  alters  the 
membrane  current  density  profiles.  For  all  three  GDLs  in  Fig.  5,  the 
current  density  near  the  land  region  becomes  greater  than  that  near 
the  channel  area,  indicating  that  the  membrane  is  not  fully  hydrated 
due  to  the  higher  rise  in  temperature.  Thus  the  PEFC  performance 
is  determined  mainly  by  the  degree  of  membrane  hydration.  Con¬ 
sequently,  in  the  three  GDL  cases,  a  lowered  cell  performance  is 
seen  in  Fig.  5  due  to  the  membrane  dehydration  effect  as  well 
as  the  electronic  contact  resistance  between  the  GDL  and  the  CL. 
The  drop  in  performance  is  31,  25,  and  18  mV  for  the  thru-plane, 
anisotropic,  and  in-plane  GDL  cases,  respectively.  The  simulation 
results  in  Figs.  4  and  5  show  that  when  both  anode  and  cathode  gas 
streams  are  fully  humidified,  the  rise  in  the  membrane  temperature 
should  be  maintained  within  2-3  °C  to  avoid  significant  membrane 
dehydration. 

The  membrane  temperature  and  current  density  distributions 
at  a  lower  current  density,  i.e.,  0.5  A  cm-2,  are  presented  in  Fig.  6 


Fig.  4.  Current  density  and  temperature  distributions  in  the  membrane  for  the  three 
GDL  cases  (/=  1.5  A  cm-2  without  the  contact  resistances  at  the  GDL/CL  and  BP/GDL 
interfaces).  Thru-plane  GDL  case:  Vcen  =  0.559  V:  anisotropic  GDL  case:  Vceu  =  0.566  V: 
in-plane  GDL  case:  Vce/I  =  0.566  V. 


where  the  anisotropic  and  in-plane  GDLs  are  considered  without 
the  effect  of  the  contact  resistance  between  the  GDL  and  the  CL.  In 
addition,  because  the  lower  current  density  (0.5  A  cm-2)  results  in 
a  lower  gas  velocity  in  the  gas  channels  under  a  constant  stoichiom¬ 
etry  situation,  a  15%  coverage  of  the  cathode  GDL  surface  by  liquid 
water  is  assumed  to  approximate  the  presence  of  water  droplets 
at  the  cathode  GDL/gas  channel  interface.  Due  to  the  low  current 
density  operation,  the  rise  in  membrane  temperature  is  not  signif¬ 
icant,  i.e.,  only  0.5  and  0.3  °C  for  the  anisotropic  and  in-plane  GDL 
cases,  respectively.  As  a  result,  the  membrane  current  density  pro¬ 
files  for  both  cases  are  similar  to  each  other,  which  indicates  that 
the  isotropic  assumption  for  the  GDL  is  reasonable  in  the  case  of 
low  current  density  operation,  e.g.,  0.5  A  cm-2. 


Fig.  5.  Current  density  and  temperature  distributions  in  the  membrane  for  the  three 
GDL  cases  (7=1.5  A  cm-2  with  the  contact  resistances  at  the  GDL/CL  and  BP/GDL 
interfaces).  Thru-plane  GDL  case:  Vcen  =  0.528  V:  anisotropic  GDL  case:  Vceu  =  0.541  V: 
in-plane  GDL  case:  Vceu  =  0.548  V. 
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Fig.  6.  Current  density  and  temperature  distributions  in  the  membrane  for  the 
anisotropic  and  in-plane  GDL  cases  (7  =  0.5  A  cm-2;  without  the  contact  resistances 
at  the  GDL/CL  and  BP/GDL  interfaces;  the  15%  interfacial  liquid  droplet  coverage 
is  assumed  at  the  cathode  GDL/GC  interface).  Anisotropic  GDL  case:  Vceu  =  0.709  V: 
in-plane  GDL  case:  Vceu  =  0.708  V. 


For  the  three  cases  presented  in  Fig.  4  (i.e.,  based  on  1.5  A  cm-2 
operation  without  contact  resistance  at  the  GDL/CL  interfaces), 
Figs.  7  and  8,  respectively,  show  the  temperature  and  liquid  sat¬ 
uration  contours  on  the  two-dimensional  cross-section  of  the  PEFC 
where  the  anode/cathode  GDLs,  anode/cathode  CLs,  and  membrane 
are  seen.  As  shown  in  Fig.  7,  while  a  considerable  rise  in  tem¬ 
perature  (4.15  °C)  occurs  in  the  thru-plane  GDL,  there  is  a  little 
difference  in  the  temperature  rises  for  the  anisotropic  and  in-plane 
GDLs,  i.e.,  2.15  and  1.35  °C,  respectively.  On  the  other  hand,  the 
effect  of  GDL  anisotropy  is  much  more  appreciable  in  terms  of 
temperature  contours.  When  the  in-plane  thermal  conductivity  is 
an  order  of  magnitude  larger  than  the  thru-plane  thermal  con¬ 
ductivity,  the  temperature  contour  in  the  anisotropic  GDL  is  more 
one-dimensional  and  varies  more  significantly  in  the  thru-plane 
direction.  Therefore,  the  isotropic  and  anisotropic  GDLs  exhibit 
vastly  different  thermal  behaviour;  this  clearly  demonstrates  that 
the  assumption  of  isotropic  GDL  thermal  conductivity  cannot  pro¬ 
vide  realistic  temperature  distributions. 


The  liquid  saturation  contours  for  the  same  three  GDL  cases  are 
presented  in  Fig.  8.  Even  with  the  same  rate  of  water  production  that 
corresponds  to  the  same  current  density,  i.e.,  1.5  A  cm-2,  it  is  clearly 
seen  that  the  liquid  saturation  profiles  differ  among  the  three  GDLs, 
being  strongly  affected  by  the  temperature  profiles  shown  in  Fig.  7. 
Due  to  the  high  rise  of  temperature  in  the  thru-plane  GDL,  most  of 
the  water  near  the  channel  regions  evaporates  and  water  accumu¬ 
lates  near  the  land  region.  On  the  other  hand,  water  accumulates 
in  a  larger  portion  of  the  GDL  for  both  the  anisotropic  and  in-plane 
GDL  cases,  due  to  their  relatively  smaller  rises  of  temperature. 

It  should  be  noted  that  two  important  phenomena  emerge  in 
the  non-isothermal  two-phase  situation  namely:  (i)  vapour-phase 
diffusion  that  is  driven  by  the  thermal  gradient  in  the  two-phase 
region  and  (ii)  the  heat-pipe  effect  due  to  the  transfer  of  latent 
heat  from  evaporation  in  the  catalyst  layer  (hot  region)  and  re¬ 
condensation  over  the  current-collecting  land  (cold  region).  Both 
phenomena  depend  on  the  thermal  gradient  in  the  GDL  between 
the  hot  and  the  cold  regions.  In  the  thru-plane  and  anisotropic  GDLs, 
the  reverse  liquid  saturation  gradient  is  observed  along  the  thru- 
plane  direction  from  the  current-collecting  land  side  of  GDL  toward 
the  CL  side  of  GDL.  This  is  due  to  the  effect  of  relatively  stronger 
vapour-phase  diffusion  that  is  driven  by  the  larger  thermal  gradi¬ 
ent.  Consequently,  the  maximum  liquid  saturation  for  these  two 
cases  occurs  near  the  surface  of  the  interface  of  the  GDL  that  faces 
the  current-collecting  land.  On  the  other  hand,  the  reverse  liquid 
saturation  gradient  is  not  seen  in  the  in-plane  GDL  because  the 
vapour-phase  diffusion  is  weak  due  to  the  small  variation  of  tem¬ 
perature  in  the  GDL.  Therefore,  the  isotropic  treatment  that  uses  the 
value  of  the  in-plane  thermal  conductivity  fails  to  predict  accurately 
the  liquid  saturation  profile  in  the  GDL  by  not  capturing  the  reversed 
liquid  saturation  gradient  that  arises  from  the  strong  vapour-phase 
diffusion. 

The  temperature  and  liquid  saturation  contours  with  reference 
to  the  CL/GDL  contact  resistance  are  shown  in  Figs.  9  and  10,  respec¬ 
tively.  From  a  comparison  of  Figs.  7  and  9,  it  is  clearly  seen  that  the 
temperature  gradient  between  the  channel  and  land  regions  in  GDL 
is  greater  in  the  case  of  Fig.  9  due  to  the  variation  of  the  contact 
resistance.  The  larger  thermal  gradient  in  Fig.  9  enhances  vapour- 
phase  diffusion  and  renders  the  level  of  liquid  saturation  near  the 
land  region  higher.  As  a  result,  the  maximum  values  of  liquid  satu¬ 
ration  in  the  thru-plane  and  anisotropic  GDLs  in  Fig.  10  are  slightly 
higher  than  the  maximum  values  in  Fig.  8.  An  exceptional  trend  is 
seen  in  the  in-plane  GDL,  where  the  rise  in  temperature  without  the 
GDL/CL  contact  resistance  in  Fig.  7  is  minimal  and  thus  facilitates 
the  condensation  of  more  water  in  the  GDL  due  to  the  lower  liquid 


(a)  Tmax :  357.3K 
Tmin  :  354.1K 

aGDL  cGDL 


(b)  Tmax :  355.3K 
Tmin  :  354.0K 

aGDL  cGDL 


(c)  Tmax :  3 54.5 K 
Tmin  :  353.9K 

aGDL  cGDL 


Fig.  7.  Contours  of  temperature  distributions  in  the  two-dimensional  cross-section  of  PEFC  where  the  anode/cathode  GDLs,  anode/cathode  CLs,  and  membrane  are  plotted 
(7=  1.5  A  cm-2  without  the  contact  resistances  at  the  GDL/CL  and  BP/GDL  interfaces). 
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(a)  Smax :  0.1853  (b)  smax :  0.1802  (c)  smax :  0.1921 
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Fig.  8.  Contours  of  liquid  saturation  distributions  in  the  two-dimensional  cross-section  of  PEFC  where  the  anode/cathode  GDLs,  anode/cathode  CLs,  and  membrane  are  plotted 
(/=  1.5  A  cm-2  without  the  contact  resistances  at  the  GDL/CL  and  BP/GDL  interfaces). 
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Fig.  9.  Contours  of  temperature  distributions  in  the  two-dimensional  cross-section  of  PEFC  where  the  anode/cathode  GDLs,  anode/cathode  CLs,  and  membrane  are  plotted 
(/=  1.5  A  cm-2  with  the  contact  resistances  at  the  GDL/CL  and  BP/GDL  interfaces). 
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Fig.  10.  Contours  of  liquid  saturation  distributions  in  the  two-dimensional  cross-section  of  PEFC  where  the  anode/cathode  GDLs,  anode/cathode  CLs,  and  membrane  are 
plotted  (/=  1.5  A  cm-2  with  the  contact  resistances  at  the  GDL/CL  and  BP/GDL  interfaces). 


268 


H.Ju  /  Journal  of  Power  Sources  191  (2009)  259-268 


saturation  pressure,  which  increases  the  level  of  liquid  saturation  in 
GDL.  Again,  there  is  considerable  difference  in  the  liquid  saturation 
profiles  between  the  isotropic  and  anisotropic  GDLs,  which  indi¬ 
cates  a  strong  need  to  consider  an  anisotropic  GDL  in  the  modelling 
and  simulation  of  PEFC. 

4.  Conclusions 

Two-phase,  non-isothermal  simulations  have  been  performed 
to  assess  the  effects  of  an  anisotropic  GDL  and  the  spatial  variation 
of  contact  resistance  between  the  CL  and  the  GDL  on  the  distribution 
of  both  temperature  and  water,  as  well  as  cell  performance.  The 
following  conclusions  can  be  reached. 

1 )  The  numerical  results  reveal  vastly  different  thermal  behaviour 
in  the  isotropic  and  anisotropic  GDLs,  and  predict  different  tem¬ 
perature  profiles  due  to  the  differences  in  the  through-plane 
and  in-plane  GDL  thermal  conductivities.  When  the  in-plane 
thermal  conductivity  is  an  order  of  magnitude  larger  than  the 
through-plane  thermal  conductivity,  the  temperature  profile  in 
the  anisotropic  GDL  is  more  one-dimensional,  and  varies  more 
significantly  in  the  thru-plane  direction.  Due  to  the  strong  inter¬ 
action  between  the  water  and  heat  transport  phenomena  in  a 
PEFC,  different  liquid  saturation  profiles  between  the  anisotropic 
and  isotropic  GDLs  are  also  captured  by  numerical  simulations. 

2)  In  terms  of  the  rise  of  temperature  in  a  PEFC,  the  maximum  rise  in 
the  anisotropic  GDL  lies  between  two  limits,  where  the  lower  and 
upper  limits  correspond  to  the  isotropic  GDLs  that  are  designed 
as  having  the  in-plane  and  through-plane  thermal  conductivities 
in  all  directions,  respectively.  Since  the  maximum  temperature 
in  the  anisotropic  GDL  leans  more  toward  the  lower  limit  that  is 
predicted  in  the  in-plane  GDL,  it  is  concluded  that  isotropic  treat¬ 
ment  that  uses  the  in-plane  GDL  thermal  conductivity  results  in 
a  smaller  deviation  in  terms  of  the  maximum  rise  in  temperature 
than  that  which  uses  the  thru-plane  value. 

3)  Both  the  temperature  and  water  distributions  are  considerably 
altered  by  the  presence  of  contact  resistance  between  the  GDL 
and  the  CL,  where  the  contact  resistances  at  the  GDL/CL  inter¬ 
faces  are  assumed  to  exponentially  increase  under  the  channel 
region  due  to  the  decreasing  clamping  pressures  that  are  applied 
under  gas  channels.  As  a  result,  the  temperature  significantly 
increases  near  the  channel  regions  and  thus  considerable  mem¬ 
brane  dehydration  occurs.  This  significantly  lowers  both  the 
membrane  current  density  near  the  channel  region  and  the  over¬ 


all  PEFC  performance.  Vapour-phase  diffusion  is  also  enhanced 
due  to  the  larger  thermal  gradient  between  the  channel  and  land 
regions  that  is  created  by  the  variation  in  the  contact  resistance 
and  results  in  a  reversed  liquid  saturation  gradient  along  the 
thru-plane  direction. 
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